Efficient linear scaling method for computing the thermal conductivity of disordered 

materials 
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An efficient order— N real-space Kubo approach is developed for the calculation of the thermal 
conductivity of complex disordered materials. The method, which is based on the Chebyshev polyno- 
mial expansion of the time evolution operator and the Lanczos tridiagonalization scheme, efficiently 
treats the propagation of phonon wave-packets in real-space and the phonon diffusion coefficients. 
The mean free paths and the thermal conductance can be determined from the diffusion coefficients. 
These quantities can be extracted simultaneously for all frequencies, which is another advantage in 
comparison with the Green's function based approaches. Additionally, multiple scattering phenom- 
ena can be followed through the time dependence of the diffusion coefficient deep into the diffusive 
regime, and the onset of weak or strong phonon localization could possibly be revealed at low 
temperatures for thermal insulators. The accuracy of our computational scheme is demonstrated 
by comparing the calculated phonon mean free paths in isotope-disordered carbon nanotubes with 
Landauer simulations and analytical results. Then, the upscalibility of the method is illustrated 
by exploring the phonon mean free paths and the thermal conductance features of edge disordered 
graphene nanoribbons having widths of ~20 nanometers and lengths as long as a micrometer, which 
are beyond the reach of other numerical techniques. It is shown that, the phonon mean free paths 
of armchair nanoribbons are smaller than those of zigzag nanoribbons for the frequency range which 
dominate the thermal conductance at low temperatures. This computational strategy is applicable 
to higher dimensional systems, as well as to a wide range of materials. 

PACS numbers: 72.80.Vp,72.15.Rn,73.22.Pr 
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I. INTRODUCTION 

Thermal conductivity of materials plays a crucial role 
in the efficiency of device applications at the nano-scale. 
In electronic applications, it is required to transfer the ex- 
cess heat effectively in order the device to work efficiently. 
For thermoelectric applications, on the other hand, a low 
thermal conductivity is essential to achieve a high per- 
formance. With the advances in the fabrication and ma- 
nipulation at the nano-scale, it has become possible to 
consider new means of thermal management like tunable 
thermal links, 2 thermal transistors 5 ' 4 and thermal logic 
gates. ' Also, it is shown that the conventional barriers 
limiting the thermal conductivity can be broken. Lying 
at the heart of a broad spectrum of applications with dif- 
ferent, if not opposite, demands, thermal management at 
the nanoscale is attracting a growing interest. 

On the other hand, carbon based materials (carbon 
nanotubes, graphene)' 1 are of special importance due 
to their very high electronic mobilities and thermal 
conduetivies. The extremely high thermal conductiv- 
ity of two-dimensional clean graphene is not suited for en- 
ergy applications. However, a further intentional damage 
(with irradiation, isotope doping,...) or a reduction of the 
dimensionality (graphene nanoribbons) of the material 



can offer the way to tune phonon conduction while pre- 
serving good electronic conductance. An ideal situation 
would be to design a material with exceptional electrical 
conduction together with a thermally insulating state. 
This turns out to be extremely challenging, and would 
require to find a way to strongly localize phonon modes, 
somehow "trapped" in a random disorder potential. The 
existence of an Anderson localization regime for acoustic 
phonons has been reported in some disordered materi- 
als 12 ' 13 , but so far not in carbon-based materials. Disor- 
dered carbon nanotubes (CNT)-based bundles exhibit a 
tendency towards a weak thermal insulating regime 14 ' 15 , 
whereas isotope disorder remains inefficient to strongly 
localize coherent phonons, even in presence of a large 
and complex underlying disorder profile. 1,1 In contrast to 
carbon nanotubes, graphene nanoribbons (GNRs) offer 
additional sources of potential phonon scattering and lo- 
calization because of the presence of irregular edges and 
enhanced chemical and structural reactivity that could 
affect low-energy phonon modes. 1. Other suggestions in- 
clude the design of heterostructures made from CNTs ls 
or pristine graphene mixed with disordered (isotope im- 
purities) GNRs 19 , or selective functionalization of GNRs 
via grafted hydrogen impurities. 20 

With these advances in nano-scale fabrication and 
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thermal management, it becomes crucial to develop new 
computational techniques which are both able to ac- 
count for the atomistic details of the systems and also 
can handle systems having sizes which are experimen- 
tally relevant. Here, an efficient linear scaling approach 
to compute coherent phonon propagation in structurally 
complex materials is described in detail. This approach 
is based on the Kubo methodology and on the exten- 
sive use of the MKRT technique, JJ which is a real-space 
(order-TV) computational framework. It has been suc- 
cessfully used for treating complex electron transport 
problems in quasi-periodic systems, two-dimensional dis- 
ordered systems in high magnetic fields 2 ", carbon nan- 
otubes 23-2 ', semiconducting nanowires 28 ' 29 or graphene- 
based materials. 30 We note that implementation of the 
Lanczos method for computing the Landauer Biittiker 
conductance of low dimensional systems has also been 
reported. 31 . 

In this paper, we first extend our recent communica- 
tion on the method 32 to a complete derivation of the 
phonon transmission coefficient starting from the origi- 
nal Kubo formula and within the framework of the har- 
monic approximation. This means that only elastic scat- 
tering (due for instance to isotope disorder) will be in- 
troduced, disregarding anharmonic effects. The study 
of the dynamical properties of phonon wave-packets will 
also be related to the thermal conductance which requires 
a phonon frequency integration over the whole spectrum. 
We then validate the method by comparison with other 
numerical approaches or analytical results, and further 
apply this method to edge disordered GNRs and discuss 
its limits. 



the Lanczos technique, one can avoid any matrix inver- 
sions, and a considerable gain in the computational effi- 
ciency is obtained, which allows the study of very large 
scale materials. 



A. Derivation of the phonon transport equations 

The vibrational Hamiltonian, taking only the harmonic 
interactions into account, is described as 



i 1 ij 



(1) 



where Ui and pi are the displacement and momentum op- 
erators for the ith atomic degree of freedom, M, is the 
corresponding mass, and <I> is the force constant tensor. 
Based on the linear response theory, the phonon conduc- 
tivity a along x direction can be obtained as 

rP r°° 

a = 9.T- 1 lim lim / dX dt e^+W* (J x {-ih\) J x (f)), 
u^o ij->-o J J 

(2) 

with being the system volume and T being the tem- 
perature. J x is the x component of the energy flux oper- 
ator J, and it can be expressed as J x = l/2Vl^2 li .{Xi — 
Xj)<f>ijUiVj, where Vj is the velocity operator and Xi 
is the equilibrium position of the atom to which the 
ith degree of freedom belongs. After neglecting terms 
like aj^aj, and d m d n , J x can be rewritten in terms of 
the phonon creation and annihilation operators as J x = 

T, m ,n J mn°L&n> where 



II. COMPUTATIONAL PHONON TRANSPORT 
METHODOLOGY 

The electronic transport theory in the linear response 
regime generally relies on the approach derived by R. 
Kubo. 33 To investigate bulk quantum phonon transport 
in disordered materials, the use of the Kubo formal- 
ism turns out to be the most natural and computation- 
ally efficient one. It has already been used for inves- 
tigating thermal transport in disordered binary alloys 
or nanocrystralline silicon. 34,35 Inspired by the MKRT 
scheme for electron transport, 21 we derive a real-space 
implementation of the Kubo formula for phonon prop- 
agation, which establish a direct computational bridge 
between phonon dynamics and the thermal conductance. 
In contrast to other implementations of the Kubo ap- 
proach, we extract the dynamical information from the 
time evolution of the wave-packet 31 ' (based on the expan- 
sion of the evolution operator on a Chebyshev polyno- 
mials basis) and simulate quantum dynamics instead of 
solving the Newtonian equations of motion. 37-39 There- 
fore, a single initial condition is enough, i.e. the initial 
atomic displacements, without any need to compute the 
time-dependent atomic velocities. Additionally, by using 



i\[X,D]\, 



(3) 



X is the diagonal matrix of equilibrium positions, D 
is the mass normalized dynamical matrix with Dij = 
$y/ y/MiMj and \n) is the nth eigenstate of D. Allen 
and Feldman have shown that Eq.(2) can be written as 



Of, 



TlT E doJm J mnJnm^m W„). 



(4) 



where fg is the Bose distribution function. Therefore, 
one has 

oo 

^ = ~ j ^^-^-Tr{[X,D]6(u-VD)[X,D]S(uj-VD)}, 
o 

(5) 

where yD = u n \n)(n\ acts as the free particle Hamil- 
tonian, and Tr{. . .} stands for the trace operator. Defin- 
ing V x = —i[X, v~D], one can write the thermal conduc- 
tance of a one-dimensional system as 



V- 



= Zc [ duj hcu^-Tr{V x 5(uj - \/D)V x 5(oj - Vd)}. 
Jo dT 

(6) 
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The thermal conductance can also be derived from the 
Landauer formalism 1 " or the nonequilibrium Green's 
function approach 41 as 



K = T^Z I dLuhuj—T(uJ) 



2tt 



(7) 



with 7~(w) being the phonon transmission function. 
Comparing these two formulas, we obtain the transmis- 
sion function as 

7» = ^-Tr{V x S(u> - \fD)V x 8{uj - VD)}. (8) 

The phonon transmission function derived here has ex- 
actly the same form as the electron transmission function 
derived from the Kubo-Greenwood formula 21 , 



Tel(E) = 



2*2 



2ir 2 h 
~1? 



■Tr{%8(E-H e i)V x 6(E-H e i)}, (9) 



where "H c i is the electronic Hamiltonian. 

We rewrite the Kubo formula in a more convenient 
form for the real-space study of the wave-packet dynam- 
ics. Expressing the 5— function as 



S(u- 



1 f°° 

d) = — toS uWB)t , (io) 

2lT J_„ 



one can then rewrite T(lv)L 2 /ir as 

dtTr [5(oj -VD)} (V x (t)V x (Q)) u 

dm- {s(u> -Vd)} (v x (t)v x (o) + v x (o)v x (tj)u) 



-,-iVDt 



and 



where V x (t) = W(t)V x U(t) with U(t) 
(. . .) u denotes the mean value of the operator over dif- 
ferent eigenstates with frequency u. The mean-square 
displacement along the x— direction over the states hav- 
ing frequency oj is written as 



X 2 (u J ,t) = ((X(t)-X(0)f) L 
where x(t) = W(t) X U(t). Hence, 



f!2) 



dt 



x '(u',t) = (V x (0)(X(0) - X(-t)) 

+(X(0)-X(-t))V x (0)) u , (13) 



and 



— x 2 (" 2 ,t) = (V x (t)V x (0) + V x (0)V x (t)) u . (14) 

Using ^x 2 (w 2 , t)\ t =o = 0, the integral in (11) can be 
evaluated as 

TH = -|> {S^ 2 D)} Hm ^{ U , t). (15) 

We consider (V x (Q)V x (t)) u = v 2 (u;)c-* /rt ', where v(u) 
and Tt r stand for the average group velocity and the av- 
erage transport time at frequency uj, respectively. With 



this approximation, the multiple scattering effects are 
taken into account in an average manner, and the mean- 
square displacement is obtained as 



X 2 (w,i) = 2v 2 r tl .t - 2i> 2 t 2 + 2w 2 r 2 .e- t / n '. 
The diffusion coefficient is defined as 



T>{u,t) 



f 



(16) 



(17) 



therefore in the ballistic regime t <C T tr and T>(ui, t) = v 2 t, 
while in the diffusive regime t 3> r tr and T>(uj,t) = 
O max (w) = 2v 2 Ttr- The transmission function in the dif- 
fusive regime reduces to 



2u>TT 



7» = — Tr {5{u 2 - D)} V max (u). 



(18) 



The phonon transport mean free path (MFP) can be ob- 
tained via 



VT tr 



2v(uj) 



(19) 



We note that, £{oS) can also be approximated by 
T(w)L/2A r c i 1 , which gives the same result as Eq.(19), be- 
cause in a quasi-one dimensional system the number of 
channels is iV c h(w) « 2w7rTr {<5(w 2 — D)\ v(u)/L. 

By computing T>(io,t), one can thus deduce T> max {uj) 
and w(w) and £(u). Noticing that 



x 2 (",t) 



Tr{(X(t)-X(0)) 2 S(u 2 -D)} 

Tr {5(lo 2 - D)} 
Tr{[X,U(t)Y5(LQ 2 -D)lX 7 U(t)}} 
Tr {6(lo 2 -D)} 



(20) 



the trace in the numerator can be calculated efficiently 
through an average over a few random phase states as 
N(i:\[X,U{t)}U(uj 2 - D)[X,U(t)]\ip), N being the num- 
ber of degrees of freedom and the bra-ket corresponds 
to a kind of projected density of states (PDOS) associ- 
ated with the vector [X, U(t)]\ip). PDOS differs from the 
normal PDOS by a factor of l/2w and is obtained by 
using the Lanczos method. Tr {S(u> 2 — D)\ is also calcu- 
lated by averaging the PDOS of through the Lanczos 
method. [X, U(t)]\ip) is calculated by using Chebyshev 
expansion of U(t), however the very low frequency com- 
ponent cannot be evaluated with a reasonable accuracy 
for large t, which will be discussed in Section II C. If 
we transform U(t) — > U(t) = c~ iDt , the corresponding 
velocity and the transport time will be transformed as 
v — > 2luv and r tr — > T tr /2a>, up to the first order approx- 
imation of the perturbation of the dynamical matrix D. 
Therefore Eq.(20) can be approximated by 



Tr{[XM(t/2uj)^S{u} 2 - D)[X,U{t/2uj)]} 
Tr{S(u 2 -D)} 



(21) 

so that we can expand IA (r) instead of U (t) in terms of 
the Chebyshev polynomials. 
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FIG. 1. (Color online) Vibrational PDOS on a single ran- 
dom phase state for a clean two-dimensional graphene system 
calculated by using the force constants given in Ref. 54. The 
Lanczos coefficients are shown in the inset. 



B. Lanczos and continued fraction methods 

Efficient computational recursion and order- N meth- 
ods have been successfully developed in solid-state 
physics starting from the pioneering work by R. 
Haydock. 42-44 The recursion methods are based on an 
eigenvalue approach of Lanczos , and rely on the com- 
putation of Green's function matrix elements by contin- 
ued fraction expansion, which can be implemented ei- 
ther in real or reciprocal spaces. These techniques are 
particularly well suited for treating disordered materials 
(alloys,..) and defect-related problems, and were success- 
fully implemented to tackle impurity-level calculations in 
semiconductors using the tight-binding approximation 4 ' 1 , 
or electronic structure investigations of amorphous semi- 
conductors, transition metals and metallic glasses based 
on the linear-mufnn-tin orbitals. 47 The vibrational DOS 
of disordered materials have also been investigated by the 
recursion method. 

Owing to the 5— function, the numerator and the 
denominator can be considered as PDOS on \ijj) and 
[X, U(t)]\ip), respectively. Given and assuming that 
[X, U(t)]\ip) is already known, we can use continued 
fraction method to calculate PDOS. We tridiagonal- 
ize the dynamical matrix D from an initial state \ip) or 
[X, U(t)]\ifi) by using the Lanczos scheme. The Lanc- 
zos coefficients a n and b n display a convenient feature: 
they converge rapidly to constants aoo and boo, respec- 
tively for a reasonable number of recursion steps. The 
PDOS is then deduced from the first diagonal element 
of retarded Green's function in the tridiagonalized repre- 
sentation. For instance, 



1 



(ip\S(uj 2 - D)\ip) = lim lm(ip\ 



1 



7T ?7->0 U! 2 + if] — D 



(22) 



ui 2 + irj — a\ — 



uj + ir} — a 2 



U) 2 + irj — aN — b 2 N Y*(u>) 



(23) 

In practice, the continued fraction is terminated after a 
certain step N by the approximation a_/v + i = and 
bN+i = boo ■ The termination can be analytically ob- 
tained as 



j 2 + irj - aoo - iyji^boo) 2 — {u 2 + 



if] - a c 



26^ 

(24) 

As an example, the normal PDOS on a single random 
phase state for a clean two-dimensional graphene system 
is given in Fig. I. It already agrees well with the DOS ob- 
tained by direct matrix diagonalization. The inset shows 
the corresponding Lanczos coefficients, a(n) and b(n). 



C. Chebyshev polynomial expansion 

The evolution of the vector [X, U(t)]\ip) is determined 
by U(t). Although U(t) has a simple form, its exact 
calculation requires the diagonalization of the dynami- 
cal matrix. Provided that all the eigenvalues and eigen- 
vectors of D are known, the matrix from of U (t) can 
be obtained by using a unitary transformation. The di- 
agonalization of a matrix requires 0(N 3 ) operations so 
it is practically impossible to handle a system with size 
N ~ I0 6 . One alternative is the Taylor expansion of 
the evolution operator. The efficiency and the accuracy 
are however difficult to fulfill at the same time with such 
an expansion, especially because a higher accuracy re- 
quires either a larger number of time steps or larger ex- 
pansion orders for a given time t. The implementation of 
the Taylor expansion for studying long time propagation 
of wave-packets in large scale disordered systems with a 
high accuracy is therefore beyond today's computational 
reach. Instead, we employ the Chebyshev polynomial ex- 
pansion approach. Formally we can expand any function 
of an operator in terms of Chebyshev polynomial-based 
operators, since the set of Chebyshev polynomials Q n 
form a complete orthogonal basis set. The expansion 
coefficients depend on the form of the expanded func- 
tion, on the time step At, and on the function domain 
[a — 26, a + 2b] which should cover the whole spectrum 
of the dynamical matrix D. Since Chebyshev polynomi- 
als are defined in the interval [—1,1], we first need to 
scale and shift the argument so that it falls within the 
range based on D' = (D — a) /2b. Finally U(At) can be 
expanded as 



and lim r j_^o(V'l( w + if] — D) can be expressed as a 



C/(At) 



-iAtVa+26-D' 



n=0 



c n (At)Q n (D') (25) 
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where 



Cn(At) 



dx'- 



1 



-iAtVa+2bx' 



7r(i + 5„,o) J_i yr^x'- 

(26) 

Thus [X,£/(Ai)]|</>) = S~ c„(At)[Jr,Qn(X>)]|V>- 
Then, the commutators are computed using the Cheby- 
shev recurrence relation bQ n +\(D') = (D — a)Q n (D') — 
bQ n - X {D') where Qo(D') = 1 and Qi(D') = (D - 
a) /2b. Recurrence relations between commutators write 
b[X, Qn+i(D')} = [A, (D - a)Q n {D')} - b[X, Q n -i{D% 
which is rewritten for convenience as 



Qn(ar') 



UD'M, \Pn) = [X, Q n (D')]\ip). (27) 



Using the well-known expression [A,BC] = [A,B]C + 
B[A, C], the commutator relationship becomes b\(3 n+ i) = 
(D - a)\p n ) - b\0n-x) + [X,D]\a n ) with \fo) = and 
= [X,D]\tp)/2b. The computation of |/3„) require the 
knowledge of vectors \a n ) = Q n (D')\ip) that will appear 
in the Chebyshev expansion of the evolution operator 
U(AT)\ip). Such \a„) vectors satisfy 



b\a„+i) = (D - a) | On) - b\a n -x), 



(28) 



with |a ) = \i>) and \a{) = (D-a)\ip)/2b. The algorithm 
thus consists of computing in parallel two recurrence re- 
lations and summing up the series 



U(MM = £ c n (At)\a n ), 



and 



n=0 



[X,U(At)]\i>) = c n(At)\f3 n 



(29) 



(30) 



In order to reach the desired accuracy, the number of 
recursion steps (iV po i y ) needs to be chosen appropriately 
depending on the evolution step and the spectral band- 
width. Once U(mAt)\tp) and [X, U(mAt)]\ip) are ob- 
tained at time mAt, then for the following evolution step 
one has 



U((m + l)At)\ip) = [X,U(At)]U(mAt)]\ip), (31) 



and 



[X,U((m + l)At)]\ip) = [X,U(At)]U(mAt)\tlj) 

+U(At)[X, U{mAt)}\tp).(32) 

Thus, the evolution of [X , U (t)]\ip) can be calculated 
step by step from any starting \ip) provided c„ in Eq.(26) 
are known. Accordingly, the algorithm scales linearly 
with the system size and the computation time. 

For U(t) = c-t^Dt. these coefficients cannot be calcu- 
lated explicitly. We therefore have to introduce a discrete 
grid and employ a numerical quadrature formula. Here, 
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FIG. 2. (Color online) Error in the modulus of the Chebyshev 
approximation for f(x) = e~ 1 ^ with x £ [0, 10 4 ], a = 5000, 
b — 2500 and ^poiy = 200. In the inset, the absolute value of 
the real parts and the imaginary parts of the corresponding 
Chebyshev coefficients are shown, respectively. 



we use the Chebyshev-Gauss grid of which the interpo- 
lating points are the N zeros of Qn(.x'): 



cos 




, fc = 0, 1, ..JV - 1 



(33) 



and the related quadrature formula is 
/(*') 



N-l 



dx - r- 

i vT 



(34) 



k=0 



The calculated expansion coefficients are plotted in the 
inset of Fig. 2 for f(x) = e"^ with x e [0,10 4 ], 
a = 5000 and b = 2500. After some steps, c n is seen 
to decay towards zero. The decay rate of the imaginary 
part is much slower than the real part, due to the fact the 
imaginary part of e~ 1 ^ is not differentiable at zero. The 
error of the approximated modulus with 7V po i v = 200 is 
shown in Fig. 2. The approximation keeps a good accu- 
racy when x is not close or equal to zero. However, the 
error around x = is very large, which indicates the time 
evolution based on expansion of U (t) is unstable. If we 
transform U(t) —J- U(t) = c~~ lDr , the Chebyshev expan- 
sion coefficients in Eq. 25 can be evaluated analytically, 



Cn(Ar) = 



dx'- 



1 



VT 



t(1 + 5„ t0 ) j_ x 

1 + n n 



(35) 



of which both the real part and the imaginary parts de- 
cay rapidly with n, resulting in a high accuracy of the 
truncation approximation for the whole spectrum. 
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FIG. 3. Time dependent diffusion coefficients T>{ui,t) for the 
system of CNT with 10.7% 14 C isotope disorder at three cho- 
sen frequencies. 
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FIG. 4. (Color online) Frequency-dependent elastic MFP (£ c ) 
in CNT with isotopic disorder (10%) obtained by using the 
Kubo method, the GF method (data taken from Ref. 16 by 
courtesy), and the analytical formula. 



III. RESULTS AND DISCUSSION 

A. Isotope-disordered CNT 

An interesting source of phonon scattering is the iso- 
tope disorder. In carbon-based materials (nanotubes, 
graphene), the controlled incorporation of 13 C impurities 
has been experimentally demonstrated. ' 5 The ques- 
tion about a possible coherent localization phenomenon 
in such disordered nano-structures has been discussed in 
Ref. 16 using the Green's function method. As a test 
case, we investigate a CNT(7,0) with 10.7% 14 C impu- 
rities to validate our numerics in comparison with the 
Green's function results 1 ' 1 and also with the analytical 
formula for the elastic MFP, 



12oiV uc iV c hM 



(36) 



where a is the length of the lattice vector in the transla- 
tional direction, iV uc is the number of atoms in each unit 
cell, p uc is the density of states per unit cell, / is the 
percentage of isotopic impurities having mass difference 
AM, and M is the average mass of the atoms (see Ref. 37 
and 51). Fig. 3 displays the evolution of the wave-packet 
dynamics for phonon modes with different frequencies. 
The linear increase of T>(uj,t) at t > observed in all 
cases indicates ballistic transport at relatively short dis- 
tances, whereas the decrease of T>(u>, t) for the mode with 
frequency ui = 900cm" 1 at t > 50 ps is a signature of lo- 
calization phenomena. The saturation of T>(uj,t) to a 
maximum value characterizes diffusive transport. From 
the saturation values, the evolution of the elastic MFP 
(£ c ~ 2£) as a function of the phonon frequency and dis- 
order features can be extracted. The results compare well 
with those obtained with Green functions (Fig. 4-main 



frame), as well as with Eq. (36) (Fig. 4-inset). Only at 
the singularities of the phonon spectrum, where disorder 
induces a broadening of states which impact on the nu- 
merical MFP, they differ more from those obtained with 
the analytical expression. 



B. Edge disordered graphene nanoribbons 

Graphene nanoribbons (GNRs) are strips of graphene 
with widths varying from a few to several tens of nanome- 
ters, depending on their fabrication processes. s In con- 
trast to the two-dimensional graphene, which is a zero- 
gap semiconductor, the narrow lateral size of GNRs en- 
tails quantum confinement effects and allows a modu- 
lation of the corresponding electronic band gap. The 
vibrational band structures are also affected by the 
confinement. 1- Two types of GNRs with highly sym- 
metric edge shapes, i.e. zigzag (ZGNR) and armchair 
(AGNR), have been predicted and experimentally ob- 
served (see Ref. 8 for a review). Here, we consider ZGNRs 
of different widths with edge disorder and evaluate the 
corresponding transport MFPs. The comparison of the 
transport MFPs and the thermal conductances of ZGNR 
with AGNR having the same width is also studied. 

We use the fourth nearest neighbor force constants for 
building the dynamical matrices. ' The ribbon widths 
are defined with the number of zigzag chains N z = 20, 40, 
and 80 for ZGNR, and the number of dimers N a = 138 
for AGNR. The relative amount of edge defects (removed 
carbon atoms at the edges) is chosen to be 10%. 

Fig. 5 (top panel) shows the frequency-dependent be- 
havior £{uS) of the zigzag ribbon with N z = 80 for 10% 
edge disorder. Large modulations of £(ui) driven by the 
underlying vibrational band structure are observed. For 
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FIG. 5. (Color online) Top panel: Transport MFP for ZGNR 
with N z = 80 (width of 17.04 nm) with disorder density 
of 10%. Bottom panel: Frequency-dependent MFP ratio 
(■n z =so/£n z =4,o (solid line) and £n z =4o/£n z =20 (dashed line). 
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FIG. 7. (Color online). Thermal conductance for the 
ZGNR(iV z = 80) (solid lines) and the AGNR(7V a = 138) 
(dashed lines) with edge disorder of 10% and various ribbon 
lengths. 



a fixed disorder strength, the MFP is found to increase 
almost linearly with the ribbon width. In Fig. 5 (bottom 
panel), the ratio /&N Z2 ar e plotted for N z = 20, 40 
and 80, keeping the width ratio the same. One observes 
that the scaling £n z =so/^n z =40 — £n z =4o/@n z =20 gener- 
ally holds. This behavior is due the fact that the scat- 
tering rate decreases with increasing width, a behavior 
previously reported for electron transport in both disor- 
dered CNTs and GNRs. 55 Since the minimum accessible 
frequency within a reasonable computation time is lim- 
ited, 2? max could not be reached for lu < 70 cm . 

Fig. 6 shows the comparison of transport MFPs 
in AGNR and ZGNR of approximately equal widths 
(17 nm) with a 10% edge disorder. The considerable 
differences in the MFPs lead to an edge shape de- 




FIG. 6. (Color online) Transport MFPs for the AGNR (with 
N a = 138) and the ZGNR (with N z = 80) at disorder density 
of 10%. 



pendent thermal conductance behavior. Different from 
the ZGNR, AGNR possesses a wide range of quasi- 
ballistic modes with MFPs as large as several fj,m around 
950 cm -1 . We note that, the ballistic and diffusive 
modes with MFPs of several hundreds of nanometers pre- 
dominate the conductance, jeopardizing any possibility 
to observe the onset of Anderson localization, as pre- 
viously discussed for small diameter disordered carbon 
nanotubes 11 ', and quantum interference effects can be ne- 
glected for samples shorter than several /im. We obtain 
the transmission according to T{uS) = N^/ (1 + L/2£(w)) 
instead of Eq. (18), in order to take the contact resis- 
tance into account. T(ui) for u> < 70 cm -1 is obtained 
by linearly interpolating between T between uj = to 
70 cm -1 , where T = 4 at u> = 0. It has been shown that 
the error caused by this interpolation is less than 1.5% 
for thermal conductance at room temperature. 

Thermal conductances for both systems are shown in 
Fig. 7. The difference between the pristine thermal con- 
ductances of AGNR (dashed lines) and ZGNR (solid 
lines) is a result of the anisotropy in the phonon dis- 
persion. Meanwhile, the phonon MFPs of AGNR are 
smaller than those of ZGNRs for the frequencies domi- 
nating the thermal conductance at low temperatures (see 
Fig. 6). These two factors cause the thermal conductance 
of AGNR to be smaller than that of the ZGNR for a fixed 
edge disorder strength and ribbon length. Although the 
difference is found to be reduced as the ribbon width in- 
creases (not shown) , coherent phonon propagation is still 
sensitive to the ribbon edge shape. 
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C. Limits of the methodology 

The phonon wave-packet dynamics is originally deter- 
mined by the time evolution operator U(t), however the 
Chebyshev expansion of U(t) has a relatively large er- 
ror for low frequencies. Therefore we employ U[t) in 
the numerical calculation. The drawback of this approx- 
imation is that, low frequency modes evolve much slower 
than high frequency modes, as a result, MFPs for those 
modes cannot be obtained within a reasonable compu- 
tation time. Also, the accuracy of this approximation 
needs to be checked when the disorder in the system is 
relatively strong. In the case of very strong disorder, the 
approximation is less accurate, and one has to use the ex- 
pansion of U(t) directly. Even though the time evolution 
based on the expansion of U(t) cannot give the correct 
information about very low frequency modes, features 
for the high frequency modes can still be extracted, pro- 
vided that the time step and the number of Chebyshev 
polynomials are chosen properly. The lowest accessible 
frequency by using U(t) can be smaller than the one from 
U (t) depending on the system. 



IV. CONCLUSION 

We have presented an efficient linear scaling approach 
to compute the coherent phonon wave-packet propaga- 
tion in real-space and to evaluate the related thermal 
conductance. The computational accuracy and efficiency 
were demonstrated for isotope disordered carbon nan- 
otubes and large width graphene nanoribbons with edge 



disorder, respectively. A strong impact of edge disorder 
profile on the thermal conductance was found, as well 
as an edge shape dependence of thermal conductance, 
opening interesting perspectives for thermoelectrical ap- 
plications. One should remark that this linear scaling 
method can be implemented without major difficulty to 
a wide range of other materials, including Boron-nitride- 
based materials or silicon-based materials (nanowires, 
super lattices, etc.). 5 ' 
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